home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 September / macformat-004.iso / Shareware City / Graphics / VideoToolbox ƒ / Utilities / Quick3 / PsychometricFit.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-07-07  |  6.1 KB  |  168 lines  |  [TEXT/KAHL]
  1. /* 
  2. PsychometricFit.c
  3. Copyright 1990 (c) Denis G. Pelli
  4. A general-purpose function that does a maximum likelihood fit of any
  5. psychometric function to psychometric data. The returned value is the level of
  6. significance at which the fit can be rejected. The degreesOfFreedom may be zero,
  7. in which case no parameters will be adjusted, but you'll get the log likelihood
  8. and significance of the supplied parameter values.
  9.  
  10. The psychometric function (which you supply as an argument) takes two arguments:
  11. a contrast and a pointer to a paramRecord. The function Weibull() is provided in
  12. Weibull.c. Others may be added as desired. It is assumed that the first
  13. parameter is the log of threshold. No assumptions are made about the other
  14. parameters, except that any that are to be iteratively fit are assumed to be of
  15. type "double".
  16.  
  17. Quick3 is a stand-alone program that uses PsychometricFit() to do the real work.
  18. I suggest you read the introductory comments at the beginning of Quick3.c
  19.  
  20. HISTORY:
  21. 4/7/90        dgp wrote it
  22. 10/29/90    dgp    tidied up the comments
  23. 11/17/92    dgp "
  24. 1/25/93 dgp removed obsolete support for THINK C 4.
  25.  
  26. SOURCES:
  27. Quick3.h
  28. LogLikelihood.c
  29. MonotonicFit.c
  30. PsychometricFit.c
  31. SortAndMergeContrasts.c
  32. Weibull.c
  33. #From Denis Pelli's VideoToolbox:
  34. VideoToolbox.h
  35. Binomial.c
  36. ChiSquare.c
  37. Normal.c
  38. SetFileInfo.c        # Used only on the Macintosh
  39. #From Numerical Recipes in C:
  40. nr.h
  41. NRUTIL.h
  42. BRENT.C
  43. F1DIM.C
  44. LINMIN.C
  45. MNBRAK.C
  46. NRUTIL.C
  47. POWELL.C
  48.  
  49. LIMITATIONS
  50.  
  51. This program uses routines from Numerical Recipes in C. They're copyrighted, 
  52. so I can't distribute them. You can order the software:
  53.     Numerical Recipes C Diskette for Macintosh $29.95
  54. and book:
  55.     Numerical Recipes in C: The Art of Scientific Computing $44.50
  56. from:
  57.     Cambridge University Press
  58.     Order Department
  59.     510 North Avenue
  60.     New Rochelle, NY 10801
  61.  
  62. Note that I have made several changes to the Numerical Recipes in C routines: 
  63. 1.In every file I changed "float" to "FLOAT", and #included nr.h. I inserted the 
  64. statement "typedef double FLOAT;" in the file nr.h. This is because the 
  65. Macintosh computes doubles much faster than floats. If you'd rather run
  66. slowly than modify your Numerical Recipes in C files, then you will need to insert 
  67. "typedef float FLOAT;"
  68. in CalibrateLuminance.c & PsychometricFit.c in order to compile those files.
  69. The rest of the VideoToolbox doesn't care.
  70. */
  71. #include "VideoToolbox.h"
  72. #include "Quick3.h"
  73. #include <nr.h>
  74. #include <nrutil.h>
  75.  
  76. #define TOLERANCE 0.001    /* fractional tolerance of log likelihood. Not critical. */
  77.  
  78. /*
  79. I hate global variables because they hide the flow of information. However,
  80. some sort of cludge is necessary to pass the extra arguments to Error(), bypassing
  81. the Numerical Recipes routines that call it, since they only pass the 
  82. parameters they know about. These static declations at least restrict the scope of
  83. these "globals" to this file.
  84. */
  85. static double Error(FLOAT *p);
  86. static dataRecord *myDataPtr;                /* for Error() */
  87. static PsychometricFunctionPtr MyPsychFun;    /* for Error() */
  88. static paramRecord myParams;                /* for Error() */
  89. static int myDegreesOfFreedom;                /* for Error() */
  90. static int iter;                            /* for Error() */
  91.  
  92. double PsychometricFit(paramRecord *paramPtr,PsychometricFunctionPtr PsychFun
  93.     ,dataRecord *dataPtr,double *logLikelihoodPtr,int degreesOfFreedom
  94.     ,double *chiSquarePtr,int *chiSquareDFPtr)
  95. {
  96.     int i,j;
  97.     FLOAT *p,**direction,ftol,fret;
  98.     dataRecord monotonicData;
  99.     double monotonicLL;
  100.     int monotonicDF;
  101.     double P;
  102.     
  103.     myDataPtr=dataPtr;        /* copy these for use by Error() */
  104.     MyPsychFun=PsychFun;
  105.     myParams=*paramPtr;
  106.     myDegreesOfFreedom=degreesOfFreedom;
  107.     
  108.     p=vector(1,degreesOfFreedom);
  109.     direction=matrix(1,degreesOfFreedom,1,degreesOfFreedom);    /* initial directions */
  110.     if(p==NULL || direction == NULL)
  111.         PrintfExit("PsychometricFit: not enough room for arrays.\007\n");
  112.     for(i=1;i<=degreesOfFreedom;i++) p[i]=((double *)paramPtr)[i-1];
  113.     for(i=1;i<=degreesOfFreedom;i++)for(j=1;j<=degreesOfFreedom;j++)direction[i][j]=0.0;
  114.     for(i=1;i<=degreesOfFreedom;i++)direction[i][i]=0.03;    /* initial step size */
  115.     ftol=TOLERANCE;    /* fractional tolerance on Error value when done */
  116.     iter=0;
  117.     
  118.     /* do it. The psychometric function is passed to Error by the global MyPsychFun */
  119.     if(degreesOfFreedom==0)fret=Error(p);
  120.     else powell(p,direction,degreesOfFreedom,ftol,&iter,&fret,&Error);
  121.  
  122.     for(i=1;i<=degreesOfFreedom;i++) ((double *)paramPtr)[i-1]=p[i];
  123.     free_matrix(direction,1,degreesOfFreedom,1,degreesOfFreedom);
  124.     free_vector(p,1,degreesOfFreedom);
  125.  
  126.     *logLikelihoodPtr=-fret;
  127.     
  128.     /* Now compute the degree of significance at which the fit can be rejected */
  129.     monotonicData= *dataPtr;
  130.     MonotonicFit(&monotonicData,&monotonicLL,&monotonicDF);    /* overwrites data with fit */
  131.     *chiSquarePtr= -2.0*(*logLikelihoodPtr-monotonicLL);    /* -2 log likelihood ratio of hypotheses */
  132.     *chiSquareDFPtr=monotonicDF-degreesOfFreedom;            /* difference in degrees of freedom */
  133.     P=PChiSquare(*chiSquarePtr,*chiSquareDFPtr);            /* significance */
  134.     return P;
  135. }
  136.  
  137. /* There is a subtlety here. I thought that I could use Powell with the whole
  138. paramRecord, yet ask Powell to only twiddle the first few parameters, figuring that
  139. even when I was asking Powell to fit only the first few parameters the other
  140. parameters would still be there in the array when the pointer to the array was
  141. passed to Error(). Alas, Powell() and its subroutines make COPIES of the array,
  142. and naturally fail to copy the non-twiddled parameters, since they don't know about
  143. them. The solution is for Error() to have its own complete copy of the paramRecord.
  144. Each time Error() is called it updates the twiddled parameters before calling
  145. LogLikelihood(), which calls the psychometric function (*MyPsychFun)().
  146. */
  147.  
  148. double Error(FLOAT *p)
  149. {
  150.     double error;
  151.     int i;
  152.     static int lastIter=0;
  153.     
  154.     for(i=1;i<=myDegreesOfFreedom;i++) ((double *)&myParams)[i-1]=p[i];
  155.     
  156.     error=-LogLikelihood(myDataPtr,&myParams,MyPsychFun);
  157.     
  158.     /* Diagnostic printout for difficult cases */
  159.     if(iter>0 && iter%50 == 0 && iter!=lastIter){
  160.         printf("Error(): Warning, %d iterations:\n",iter);
  161.         printf("logAlpha %5.2f, beta %5.1f, gamma %5.2f, delta %6.3f -log likelihood %9.0g\n"
  162.             ,myParams.logAlpha,myParams.beta,myParams.gamma,myParams.delta,error);
  163.         lastIter=iter;
  164.     }
  165.     return error;
  166. }
  167.  
  168.